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ABSTRACT 

We present a general linear dispersion relation which describes the coupled 
behavior of magnetorotational, photon bubble, and convective instabilities in 
weakly magnetized, differentially rotating accretion disks. We presume the 
accretion disks to be geometrically thin and supported vertically by radiation 
pressure. We fully incorporate the effects of a nonzero radiative diffusion length 
on the linear modes. In an equilibrium with purely vertical magnetic field, the 
vertical magnetorotational modes are completely unaffected by compressibility, 
stratification, and radiative diffusion. However, in the presence of azimuthal 
fields, which are expected in differentially rotating flows, the growth rate of 
all magnetorotational modes can be reduced substantially below the orbital 
frequency. This occurs if diffusion destroys radiation sound waves on the 
length scale of the instability, and the magnetic energy density of the azimuthal 
component exceeds the non-radiative thermal energy density. While sluggish 
in this case, the magnetorotational instability still persists and will still tap 
the free energy of the differential rotation. Photon bubble instabilities are 
generically present in radiation pressure dominated flows where diffusion is 
present. We show that their growth rates are limited to a maximum value which 
is reached at short wavelengths where the modes may be viewed as unstable slow 
magnetosonic waves. We also find that vertical radiation pressure destabilizes 
upward propagating fast waves, and that Alfven waves can be unstable. These 
instabilities typically have smaller growth rates than the photon bubble/slow 
modes. We discuss how all these modes behave in various regimes of interest, 
and speculate how they may affect the dynamics of real accretion disk flows. 
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1. Introduction 

The physical state of the radiation pressure dominated, innermost regions of accretion 
disks around black holes has been uncertain ever since the early days of accretion disk 
theory. Standard alpha disk models in which the viscous stress is assumed to scale 
with the radiation pressure are subject to thermal and viscous instabilities (Lightman & 
Eardley 1974, Shakura & Sunyaev 1976). These instabilities are sensitive to the assumed 
prescription for the anomalous viscosity (e.g. Piran 1978), so it remains unclear whether 
or not they are actually present in real flows. If they are, the accretion flow may adopt a 
radically different state from that usually envisaged in thin accretion disk theory, such as 
the multiphase equilibrium recently proposed by Krolik (1998). 

In addition to these secular instabilities, dynamical instabilities also exist. First and 
foremost, a differentially rotating flow with a negative angular velocity gradient and an 
initially weak magnetic field is unstable to the magnetorotational instability (MRI, Balbus & 
Hawley 1991). The turbulence resulting from this instability is currently the most plausible 
mechanism known for generating the anomalous viscosity required in accretion disk models. 
Dissipation of magnetohydrodynamic waves excited by this turbulence by photon diffusion 
and photon viscosity has recently been examined by Agol & Krolik (1998). Photon diffusion 
might also affect the linear development of the instability itself in a laminar, radiation 
pressure dominated flow, but this issue has not been examined previously. 

Gammie (1998) has suggested that the overstable photon bubble modes discussed by 
Arons (1992) in the context of X-ray pulsars also exist in radiation dominated accretion flows 
in general. However, his instability analysis was limited to studying a static equilibrium 
where the effects of differential rotation were entirely neglected. Pietrini & Krolik (2000) 
have recently investigated convective instabilities in unmagnetized, radiation pressure 
dominated, differentially rotating flows. They assumed a constant vertical density profile, 
which necessarily leads to an unstable entropy gradient. Whether or not such unstable 
gradients exist in reality depends on the vertical dissipation profile, which in turn depends 
on the structure of the MHD turbulence. 

These dynamical instabilities may all play a role at some level in the dynamics and 
thermodynamics of the radiation pressure dominated portion of accretion disks. It is the 
purpose of this paper to provide a unified description of all three by deriving a general linear 
dispersion relation which incorporates them all. Such an analysis will hopefully prove to be 
a useful guide to numerical simulations which explore the nonlinear development of these 
instabilities and their effects on the resulting turbulent state of the inner accretion disk. 

This paper is organized as follows. In section 2 we discuss our basic equations and 
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assumptions. Then in section 3 we focus exclusively on the MRI by deriving its dispersion 
relation in the absence of vertical stratification so that the other instabilities are suppressed. 
We are therefore able to study just the effects of photon diffusion on what is probably 
the most important of these instabilities. In section 4 we then present the full dispersion 
relation for a vertically stratified, differentially rotating medium, and discuss its solutions. 
We discuss the relevance of these solutions to astrophysical accretion disks in section 5 and 
then summarize our conclusions in section 6. 




2. Equations and Assumptions 

In Eulerian coordinates in an inertial frame, the equations of radiation MHD which we 
use in this paper are 

^ + V-(pv) = 0, (1) 
VvJ = -Vp-pV$ + -^(V x B)xB + ^F, (2) 

dE 4 

— + v VE + -EV-\ = -V-F, (4) 

ot 3 

F = --^VE, (5) 
Vx(vxB). (6) 



and 



dB 

~dt 

Here p is the fluid mass density, v is the velocity, p is the gas pressure, B is the magnetic 
field, 7 is the gas adiabatic index, E is the radiation energy density, F is the radiative flux, 
and c is the speed of light. We have followed the notation of Stone, Mihalas, & Norman 
(1992) with regard to variables associated with the radiation field. In order to facilitate 
comparison with the work of Gammie (1998), we note that his p g is our p, his J is cE/4tt, 
and his H is F/47T. 

We neglect the self-gravity of the fluid, and assume that the gravitational potential $ 
originates from the central mass. We completely neglect relativistic effects. We also assume 
a pure electron scattering opacity n es , neglecting absorption opacity entirely. This should 
be a good first approximation to the inner parts of accretion disks around black holes 
which are electron scattering dominated. Note, however, that this has the consequence that 
the gas and radiation are coupled together purely by momentum exchange. They do not 
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exchange heat, which is why the gas energy equation (|) simply reduces to the condition 
for adiabatic flow. This approximation still recovers the basic MRI and photon bubble 
instabilities, and so we believe captures the physics that we wish to explore here. 

We have closed the radiation moment equations by assuming that the radiation field 
is close to isotropic, so that the stress tensor is diagonal with elements given by one third 
the radiation energy density. We are therefore neglecting the effects of photon viscosity 
considered for example by Agol & Krolik (1998). The radiation momentum equation (|5|) 
deserves special comment. We have assumed a simple diffusion form, and have in particular 
neglected time derivative and velocity-dependent terms that, e.g. Stone et al. (1992) and 
Gammie (1998) have retained. The full equation is 

_(_ + v . vF + Fv . v j = __ v£ _^ F . (7 ) 

We are interested in perturbations with wavelengths A smaller than the vertical disk scale 
height H, and angular frequencies u or growth rates which are of order the local angular 
velocity of the disk f2 (MRI, Balbus & Hawley 1991) or possibly a factor (2nH/ A) 1//2 larger 
(photon bubbles, Gammie 1998). All the terms on the left hand side of equation (|7|) will be 
negligible provided 

where r H = n cs pH is the scattering optical depth across a scale height, and c s is the total 
(gas plus radiation pressure) sound speed. 

Well outside the innermost stable circular orbit, the scattering depth in the inner, 
radiation pressure dominated zone of a standard a-viscosity disk (Shakura & Sunyaev 
1973)[] is given by 

3/2 

(9) 



_! / l yWRX 

t h ~ 1 x a rj\ — 

VL Edd / \r g J 



Here rj ~ 0.1 is the radiative efficiency, L/L Edd is the luminosity scaled with Eddington, 
and r g is the gravitational radius of the central object. The scattering depth is therefore 
always significantly greater than unity except perhaps in the innermost parts of disks 
accreting near the Eddington limit. Radiative equilibrium, hydrostatic equilibrium, and the 



J It is almost certainly the turbulence generated by the MRI which is responsible for "viscous" stresses 
in accretion disks, and therefore a value of a. We are nevertheless using a to give us some handle on the 
ambient conditions in a fictitious, laminar accretion disk prior to the development of any instability. The 
reader should therefore view these scalings with appropriate caution. 
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assumption that the viscous stress is a times the total pressure leads to a ratio of the speed 
of light to the sound speed c s ~ (4i?/9p) 1//2 given by 



c 

— ~ (XTu- (10) 

C s 



Hence conditions (§) are almost always satisfied in the interior of accretion disks, provided 
we do not consider wavelengths that are extremely small. The only exception is for 
near-Eddington disks in the innermost few gravitational radii, where our Newtonian 
treatment of the magnetohydrodynamics breaks down anyway. We therefore choose to 
neglect the left hand side of equation (0), which therefore reduces to the diffusion equation 
(H). The fact that Gammie (1998) found flux limited electromagnetic waves as solutions of 
his dispersion relation is a reflection of the fact that he retained these terms. In most cases, 
these terms have negligible effect on the photon bubble modes, which are at much lower 
frequency.0 Note that equation fl5|) allows us to write the radiation pressure very simply 
in terms of a gradient in the radiation energy density E, so the momentum equation (0) 
becomes 

P (?r + v " Vv I = -V (p+ -e] - pV$ + -?-(V x B)xB. (11) 



s dt J V 3 / 4vr 

Eliminating F using equation (^j), we may write the radiation energy equation (|j) as 

dE 4 ( c \ 

+ v . VE + -EV • v = V- VE . (12) 

Ot 3 ) 

In contrast to the radiation momentum equation (^), the left hand side of this equation can 
be neglected only if 

lu 4tt 2 f c\ Atx 2 fH\ 2 , v 2n fH\ , . 

a 5^ k) ~ ~ a It) and ~c <K ir H (t) • (13) 

Because a may be quite small and t# is not very small, these conditions may not be 
satisfied, so we must retain the left hand side of the radiation energy equation (||). In fact, 
it is neglecting the right hand side of this equation (photon diffusion) which returns us to 
the standard ideal MHD equations with gas and radiation pressure. When the conditions 
(|j~3l) start to become satisfied, photon diffusion starts to modify standard MHD. It is this 
effect which gives rise to photon bubbles, and which can modify the MRI. We therefore 
choose to retain the left hand side of equation (Q), although the right hand side can be 
dominant. 



2 Indeed, we have numerically recomputed all the results in this paper with the terms on the left hand 
side of equation (Q) included and found no significant differences for realistic physical parameters. 
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2.1. Equilibrium 

We consider a narrow annulus of the disk which, in equilibrium, is stationary and 
axisymmetric with a purely azimuthal flow velocity v = RQ<p. (We assume cylindrical polar 
coordinates {R, </>, z} throughout this paper.) Following the usual analysis of the MRI as 
well as previous analyses of photon bubbles, we assume the existence of a locally uniform 
magnetic field. Because we are considering a differentially rotating system, we assume this 
field has no radial component, which would otherwise be sheared into a time-dependent 
azimuthal field. Hence B = B^<j> + B z z. With the exception of the external gravitational 
potential $, we assume that all radial gradients are much smaller than vertical gradients of 
the same quantities, i.e. that the disk is geometrically thin. The radiative diffusion equation 
(||) then implies that the equilibrium radiative flux is mainly in the vertical direction, and 
thus the radial component may be neglected. We will allow for the possibility of radial fluid 
gradients creating small departures from Keplerian flow by retaining a general form for the 
angular velocity profile Q(R) in our dispersion relation. 

With these assumptions, the only nontrivial equilibrium equations are those expressing 
vertical hydrostatic balance, 

_I* + «= F ._, = , (14) 
p oz c 

radiative equilibrium, 

BF 

and radiation diffusion, 

F, = (16) 

6K cs p OZ 

Here g = d^/dz is the magnitude of the local vertical gravitational acceleration produced 
by the tidal field of the central mass. 

This equilibrium is characterized by a number of parameters which are important in 
describing wave modes. First are the gas and radiation sound speeds, defined as 

c g = (2?) and c r = ^ , (17) 

respectively. Then there are the Alfven speeds associated with each component of the 
magnetic field, 

^ = (l^F and VAz = WpJt*- (18) 
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Internal hydro dynamic modes are driven by equilibrium gradients with associated oscillation 
frequencies. The square of the epicyclic frequency is given by 



K 



J__d_ 



(19) 



It is convenient to define the square of the Brunt- Vaisala frequency separately in the gas 
and radiation as 



N 2 



9 



pc 2 g 



dp 

dz 



and 



N' 



1 fdE' 



dlnp 

dz 

dlnp 



3pq? \dz J dz 
respectively. Then the square of the total Brunt- Vaisala frequency is 



N 2 



9 



d 



p(cl + q?) dz 



P + - 3 E 



dlnp 

dz 



clNl + c 2 N, 



cl + c 2 



(20) 
(21) 

(22) 



Finally, there is a characteristic length scale associated with radiative diffusion, given 
by how far photons diffuse on an orbital time. This can be used to define a diffusion 
wavenumber as 

'?>KK cs p\ 1/2 



k, 



diflf 



(23) 



2.2. Perturbation Equations 

Linearizing our equations of motion about the equilibrium, and assuming a WKB 
space-time dependence oc exp[i(J k R dR + J k z dz — out)}, we arrive at the following 



perturbation equations: 

-iu>5p + ip\t-5v + 5v z ^- = 0, (24) 

dz 

- iup5v R - 2n P 5v^ = -ik R (dp + ^5E^j - i^B^B* + B Z 5B Z \ + i^-B z 6B R , (25) 

- iujpSvtp + ^pSvn = i-^BJB^, (26) 

ZS I 47T 

/ 1 \ k 
- iup5v z = —ik z \ 5p + - g5p - i-^B^SB^, (27) 

dp 

— icu5p + Sv z — = —i'jpk ■ 8v, (28) 

dz 

dE 4 

-iu8E + 5v z —+i-Ek-5v = -ik-5F, (29) 
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zc 

5F R = -- k R 8E, (30) 

<5K es p 

SF^ = 0, (31) 

5F Z = -F z 6 -P - -^k z 5E, (32) 
p 3K es p 

- iuj5B R = ik z B z 8v R , (33) 

- iuSBf - R -^ B R = -iB^k ■ 5v + ik z B z 5v,p, (34) 

and 

— iuj5B z = —ik R B z 5v R . (35) 

The radial and vertical induction equations together guarantee that the divergence of the 
magnetic field perturbation vanish. Hence one or the other can be replaced with 

ik R 5B R + ik z 5B z = 0. (36) 

Equations (124]) and (|2~8"D- (j32"D can be solved to yield a simple expression for the 
perturbed total pressure, viz. 

dp + -8E = [c + c r )5p - - c N + c v N r , (37) 

3 V s Wdifr / 9 V B s uJdis J v 3^ diff p 



where 



co> diff = uj + (38) 



k 2 

^diflf 



Each term on the right hand side of equation (|37| ) has a simple interpretation. The first is 
the acoustic pressure perturbation in response to rarefaction or compression. The second 
is the buoyancy force (note that iSv z /u is simply the vertical Lagrangian displacement of 
a fluid element). The last term is the radiation pressure perturbation arising from the 
enhanced (diminished) optical depth in compression (rarefaction) regions. The middle term 
can give rise to convection, while the last term is responsible for photon bubble modes. 



The twelve equations (p4|)-(|34|) and (j36|) relate twelve perturbation variables. Eight of 
these equations have an explicit uj dependence, and so we therefore expect an eighth order 
dispersion relation. However, as we show in Appendix A, one of these modes is in fact a 
zero frequency solution which merely represents a perturbation to the equilibrium. The 
interesting physics therefore lies in a seventh order dispersion relation. In section 4, we will 
present this full dispersion relation, but we will first discuss in section 3 the simpler case of 
the MRI in a non-stratified equilibrium. 
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3. The Magnetorotational Instability in the Absence of Stratification 

In the fictitious laminar flow we are considering, it is presumably the MRI which 
generates the turbulence whose dissipation, combined with cooling, sets up the thermal 
profile which then might drive photon bubble and convective modes. We therefore 
begin by considering the MRI in isolation, but in a flow which is still hot enough that 
radiation pressure dominates gas pressure. In the absence of equilibrium stratification, the 
perturbation equations can be solved to yield the following dispersion relation: 

D ms D BU + #AA+ - ^) - | ^ 4 = 0. (39) 

Here 

^ = ^ - k y Az , (40) 

Ans = kl U + ^ + V\\ (41) 

V s ^difr v 

are magnetosonic dispersion terms, and 

£bh = ^ - ^ - 4tfk 2 z v 2 Az (42) 

is the MRI dispersion relation (if set to zero) in a non-stratified, incompressible medium 
(Balbus & Hawley 1991). 

The dispersion relation (|39| ) describes how the MRI is coupled to magnetosonic modes 
when allowance is made for compressibility and radiation. Without radiation, it reduces 
to equation (64) of Blaes & Balbus (1994) for = 0, which is all they considered. The 
conclusion of those authors remains valid here: provided the equilibrium azimuthal magnetic 
field is sufficiently subthermal, with respect to either the gas or radiation energy densities, 
then the MRI remains essentially unchanged. Mathematically, this is because the Alfven 
speed is then much smaller than the relevant sound speed. For the MRI, \u\ ~ f2 ~ kv^ z , 
and the dispersion relation reduces to 

U + —A D BH = 0, (43) 
V g ^difr / 

so that the MRI mode frequencies remain unchanged from their incompressible form. 
Physically, in subthermal magnetic fields, the relevant wave speed for the MRI (the Alfven 
speed) is much slower than the sound speed, so the associated fluid motions are very nearly 
incompressible. 
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The question of how subthermal the azimuthal field must be in order not to affect the 
MRI is determined by radiative diffusion. From the form of the dispersion relation (|3^), it 
is clear that the MRI will be unaffected provided 



A sufficient condition for this to be true is that VA<t> <C c g . However, in a radiation pressure 
dominated region of an accretion disk, it is reasonable to consider azimuthal fields for which 
c g < t>A</> < c r . In this case the MRI can be modified if the characteristic wavenumber for 
instability £1/va z exceeds the diffusion wavenumber kdig, i.e. photons have time to diffuse 
across a wavelength in an orbital time. This reduces the effective sound speed, making 
the fluid more compressible. The criterion that the MRI be unaffected is more stringent: 
va<p < c r (k diS v Az /tt). 

Compressibility acts to reduce the growth rate of the MRI in the presence of azimuthal 
fields, but it does not remove the instability (Blaes & Balbus 1994). This is illustrated in 
Fig. 1, which depicts the maximum growth rate as a function of azimuthal field strength 
and diffusion wavenumber. We have restricted consideration to vertical (k^ = 0) modes, as 
these are expected to have the fastest growth rates, and we consider an equilibrium with 
negligible gas pressure (c g = 0). At high values of kdisVAz/Q, photon diffusion is negligible, 
and we recover the MRI behavior for an ideal compressible fluid. Subthermal azimuthal 
fields then produce no effect on the growth rate of the instability, but when the field 
starts to become thermal in strength, there is a modest reduction in the growth rate. For 
kdifiVAz/Q ~ 1, photon diffusion is important. In this regime, azimuthal fields can reduce 
the growth rate quite substantially if VA<f> ^ c r (k,imv Az/^)- For the parameters chosen in 
the figure, this corresponds to va^/vaz <^ 1 for fcdifr^Az/^ = 0.1. Higher values of kdisVAz/Q 
(less diffusion) require higher azimuthal field strengths to affect the instability. 

We can estimate how small the MRI growth rate becomes when diffusion is important 
by setting k^ = c r = in equation (|39|) and then solving it in the limit c g — > 0. The MRI 
then satisfies 



3 We do not have a rigorous proof that the actual stability criterion remains unchanged in our case. 
However, we suspect that this is true based on our numerical calculations illustrated in Fig. 1. In addition, 
the dispersion relation ( |38| ) admits a zero frequency root for ui when k 2 v Az = —R(d£L 2 /dR) (see section 
4) , which is the condition for the marginally stable wavenumber of the MRI in incompressible MHD in the 
absence of stratification. 




(44) 




(45) 



viz (flg + *M*)+*i*(rf + *M.)' 
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For 3> vaz, the unstable growth rate is therefore reduced to ~ fi(c g /t>A</,). 

Note that regardless of photon diffusion, the MRI is virtually unaffected by 
compressibility provided the field is sufficiently vertical, i.e. Va</> <C v\ z (e.g. the 
Vh<t>/vAz = 0.1 curve shown in Fig. 1). This is because a flow with a purely vertical 
equilibrium field will have vertical (/cr = 0) MRI modes identical to the incompressible 
case, as the dispersion relation (|39"D permits -Dbh = as a solution. 

To summarize, the MRI is essentially unaffected by photon diffusion, unless the 
azimuthal component of the field is comparable to or larger than the vertical component, 
and the vertical component of the field is weak enough that photon diffusion is important 
on the scale of an unstable wavelength ~ v\ z /tt. While the instability still persists, even 
subthermal azimuthal field components can then reduce the growth rate to ~ ^{cg/v^) if 
condition (f44|) is violated. 

3.1. The Relation Between the MRI, Alfven and Slow Waves 

In the next section we will examine stratified equilibria, which turn out to be subject 
to a number of unstable modes which can couple to each other in complicated ways. To 
help disentangle this complexity, it is helpful to briefly consider the relation between the 
MRI and classical MHD waves. 

Normally, in the absence of radiative diffusion, the MRI can be viewed as a differential 
rotation driven destabilization of the slow magnetosonic mode in a weakly magnetized 
medium (Balbus & Hawley 1998). The distinction between slow and Alfven modes is 
delicate, however. In the limit where the magnetic field becomes very weak, the slow 
magnetosonic mode becomes incompressible and has a dispersion relation which is the same 
as that of Alfven waves. It is in fact identical in this limit to the "pseudo- Alfven" mode 
which exists in incompressible fluids. The prefix "pseudo" is used because, although the 
mode is degenerate with the Alfven mode, the eigenfunction is still different. Alfven waves 
have velocity perturbations which are perpendicular to the plane of the wave vector and 
equilibrium magnetic field. Slow, or pseudo- Alfven, waves have velocity perturbations in 
the plane of the wave vector and equilibrium magnetic field. For example, if = in our 
geometry, Alfven waves would involve only 5v^, whereas slow waves would have nonzero 
5vr and 5v z . 

When the Coriolis forces associated with rotation are included, radial and azimuthal 
motions are coupled together, so the distinction between slow and Alfven modes is blurred 
further. It is therefore perhaps accidental that the stable modes to which the MRI converts 
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are normally the slow modes at high wavenumber. Figure 2(a) illustrates this fact in the 
case of negligible radiative diffusion. However, when radiative diffusion is important, the 
MRI can convert to the Alfven wave [Fig. 2(b)] or the slow wave [Fig. 2(c)]. 

This would only be of academic interest were it not for the fact that it is the slow mode 
which continuously transforms at low wavenumbers into the unstable photon bubble mode 
in static, stratified equilibria (Gammie 1998). The relationship between this instability and 
the MRI in linear theory therefore depends on whether the Alfven or slow mode is more 
closely tied to the MRI. Also, as we will show in the next section, the Alfven mode is itself 
unstable in stratified equilibria. 



4. The Dispersion Relation in the Presence of Vertical Stratification 



With considerable algebra, the general perturbation equations from section 2 can be 
used to derive the following seventh order dispersion relation: 



D ms D BH + klv\ z v 



D 



BH 



z u Az u A<j> 

kl ) 



klR 



dQ 2 
~dR 



/ ki 



k*k\F z + f_ i 



klF z 



+i2ttujk z k R VAzVAct ) 



k z k 2 F z 



\ldp 



3wdifrP/ pdz 



— i 



k 2 z F z I dp 



3u diS p 3uj diS p p dz 



~ ~ 2ik z g 



+k 2 u 2 ( clNl + 



CO 



^diflf 



v M , k 2 R k z u 2 F z uj 2 gdp 2 2 2 



3u; diff p 



p dz 



0. 



(46) 



The first three terms are identical to the non-stratified dispersion relation ( p9|) discussed 
in the last section. All other terms arise from the stratification. Note that in order for 
the WKB approximation to be consistent, we must have \k z \ 3> \dlnp/dz\, so some of the 
terms with explicit dependence on the vertical density gradient should probably be taken 
to be negligible. In a radiation pressure dominated accretion disk, the WKB approximation 
also demands that 

" HrJ '~ (47) 



k > 4 - 
cf 



cc; 



The dispersion relation (TO) reduces to those of previous authors under various limits. 
First, in the incompressible limit where (c^+LX^/u&g) — > oo, we recover the MRI dispersion 
relation in the presence of vertical stratification (Balbus & Hawley 1991), albeit slightly 
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modified by radiative diffusion: 

/^7\r2_L,.,„2 at? /,.,..„M 

u 2 - AQ 2 k 2 z v 2 Az = 0. (48) 



k 2 ~ 4 



fcj I c? + a;c2/a; d i ff 



g 

In addition, if we first assume that gas pressure is negligible (c 2 —> 0), adopt the ansatz that 
uj ~ k l l 2 and then take the short wavelength limit k — > oo, we obtain two roots given by 

= -*(¥)(!)'■ < 49 > 

This is identical to the asymptotic photon bubble instability dispersion relation obtained 
by Gammie (1998, eq. 41), adapted to our geometry. In spite of the presence of rotation, 
the photon bubble instability survives unscathed in this limit. 

Finally, if we set the equilibrium magnetic field equal to zero, the equilibrium density 
to be constant, and the gas pressure to be negligible, then after cancelling two resulting 
zero modes, the dispersion relation reduces to fifth order: 

o <! • k 2 c 4 / 2 Ak 2 E\ „ .k 2 cn 2 2 (Ak 2 n 2 E ol2 2 \ k z klgcK 2 n . . 
3cu 5 + i to 4 - 3k 2 + tu 3 -i tu 2 + — 3klg 2 uj + ^^ = q. (50) 

K es p \ 3p J K es p \ 3p J K es p 

This equation resembles the dispersion relation presented by Pietrini & Krolik (2000) when 
the dissipation terms that they included are neglected. Their dispersion relation still has 
several additional terms compared to equation (|50"D , however. Each of these terms are small 
and arise from the left hand side of the full radiation momentum equation (|7|), which we 
neglected. 

Despite the complexity of the general dispersion relation, the MRI turns out to be a 
robust instability. An analytic indication that this is so can be found by just setting uj = 
in equation (f46f) , which leads to the following equation: 



**■ K + fl 3 [ ii w a + ( 9 - k W) -M - k '4 - ***** = °- 

(51) 

If the gas Brunt- Vaisala frequency N s — 0, this equation shows that there will exist a zero 
frequency root if k < k crit = [—R^d^/dKjY^/vAz, the usual critical wavenumber for the 
onset of the MRI. Even if N g ^ 0, there will be a zero frequency root for kn = and 
k 2 v\ z = —R{dVt 2 /dR). While not a rigorous proof, these facts suggest that the critical 
wavenumber for the onset of the MRI is unchanged by the effects of radiation. Hence 
the instability criterion for the MRI is still that the angular velocity decreases outward 
in magnitude. These expectations are born out by our numerical calculations. Radiative 
diffusion does not affect the stability criterion, although it can dramatically affect the 
growth rates, as discussed in section 3. 
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4.1. Photon Bubble Modes and Other MHD Wave Instabilities 



w + : J " k l< z )[^ ~ " 2 k 2 (4 + < z + v 2 A(j> ) + k 2 z v 2 Az k 2 c 2 g ] = 0. (52) 



In the high wavenumber limit k — > oo, the dispersion relation can be factored to give 
k 2 c 

The first factor is a strongly damped radiation diffusion mode, the second describes Alfven 
waves, and the third corresponds to the fast and slow magnetosonic waves. Rotation is 
irrelevant at these high wavenumbers. We have lost the radiation effects on the fast and 
slow modes because in this limit the wavelength is much less than the diffusion length scale. 

It is worthwhile examining the behavior of the MHD wave modes to higher order in 
fc -1 , because they all turn out to be unstable for sufficiently high vertical photon fluxes. 
If, as we have been generally assuming, the radiation and magnetic energy densities in the 
plasma are much larger than the gas thermal energy density, then the fast and slow modes 
have frequencies given by 



u ~ ±k{v 2 Az + v 



2c« + v 2 



J A<t>) 

and 



Cr I T?^L + VU I =F 



k R 2 2 \ k z k R F z v Az 



k 2 Az ^ A V ^p(v 2 Az + v 2 )^ 



(53) 



— (1)^(1). 

respectively. The lowest order correction to the Alfven mode is of order k' 1 , and this 
correction is real, implying that any instability or damping is at still higher order. As we 
shall see, instabilities in the Alfven waves peak at finite wavenumbers. 

The first imaginary term in equation (^) represents damping of the fast modes by 
radiative diffusion. This term agrees with the optically thin damping rate derived by Agol 
& Krolik (1998), although the waves need not be optically thin for this rate to apply. The 



last term in equation (|53|) shows that a sufficiently large vertical radiative flux can overcome 
the damping and destabilize the fast waves. The instability criterion can be obtained by 
comparing the last two terms of equation fl5"3"D: 



F - >3 " c '( i+ St)^ (i ' L+i ' i * )1/2 - (55) 

Crudely, if the speed at which photons are diffusing upward in the equilibrium exceeds the 
Alfven speed, then fast modes propagating upward at some nonzero angle to the vertical 
will be unstable. For large photon fluxes, the growth rate of this instability is 

^>^p#)(§) 2 m 
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The instability criterion is most easily satisfied, and the largest growth rates are achieved, 
for vertical magnetic fields. Azimuthal fields tend to be stabilizing. 



Equation (p4|) reveals that slow modes propagating downward at some nonzero angle 
to the vertical are also unstable at high photon fluxes. These modes have the fastest 
growth rates, and at lower values of k become the photon bubble modes described by 
Gammie's (1998) dispersion relation (f49|). Adopting Gammie's (1998) photon bubble ansatz 
of to ~ k 1 ! 2 and retaining only the leading order terms as k — > oo in the general dispersion 
relation, we find 

2 [,2 2 .k z k^ es F z . 3K es pc r k z ^ ( B z 



^=[K4-i-^^-i.-^ u-) . (57) 



The middle term gives the photon bubble dispersion relation ( [49]) when use is made of the 
hydrostatic equilibrium condition ( |i4|) in the radiation dominated limit. The other terms 
compete with the photon bubble term at low and high wavenumbers, leading us to expect 
that the photon bubble dispersion relation (|49|) will only hold in the range of wavenumbers 
given by 

Beasts 4. (58) 

c 2 g c 2 

A condition equivalent to this was also noted by Gammie (1998). When the last inequality 
breaks down at high wavenumbers, the slow mode dispersion relation is better described by 
equation (|54l) . 



Figure 3 presents a numerical solution to the full dispersion relation (pE6| ) which 
illustrates these high wavenumber instabilities and compares their growth rates to the 
analytic expressions above. In addition to the fast and slow modes, the Alfven modes also 
appear to go unstable for large photon fluxes. In contrast to the fast and slow modes, 
the Alfven wave growth rates tend to zero as k — > 00. Rotation directly affects the 
unstable growth rates of the Alfven modes, and also affects the fast and slow modes at low 
wavenumbers. We describe this, and the coupling to the MRI, in the next subsections. 



4.2. Zero Azimuthal Field Case 



We now return to low wavenumbers where rotation is important. We first consider the 
case where the equilibrium field is vertical (B^ = 0). The general dispersion relation ((46|) 
then reduces to 
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+k^ ( c\Nl + -^c*nA + k * k f F ^ = o. (59) 

Note that -Dbh = is a solution of this equation when kn = 0, so we immediately conclude 
that in the absence of an azimuthal field, the vertical MRI modes are completely unaffected 
by compressibility, stratification, and radiation. This is not surprising, as the vertical MRI 
modes involve horizontal motions only, and therefore do not involve any excess forces from 
the vertical radiative flux. 

Equation ( j59|) is complicated, with many parameters. For numerical computations, we 
restrict consideration to cases where the gas pressure is negligible, so that c g and N g are 
both zero. In addition, we neglect convection and couplings to gravity modes by considering 
only isentropic equilibria, so that iV r = 0. Under these assumptions, we may use equation 
(|I4|) to write the acceleration due to gravity and the equilibrium vertical density gradient 
in terms of the photon flux, viz. 

9 = ^F z and ^ = ~^F Z . (60) 

Hence the only new parameter is F z . 

Figure 4 presents numerical calculations of the unstable mode growth rates in the 
diffusive regime at low wavenumber, for various values of photon flux. The introduction 
of stratification causes the slow (photon bubble) and Alfven wave instabilities to appear 
and have increasing growth rate as F z gets larger. The fast wave instability appears only 
after a critical threshold in F z is reached [cf. Fig. 4(b)]. This is quantitatively consistent 
with our instability criterion (1551), which predicts that the fast wave will be unstable for 



F z /{2>pv\ z ) > 125 for the parameters chosen in Fig. 4. Fig. 4(c) should be compared with 
Fig. 3 which has similar radiation parameters. Note that for the value of k^m chosen in 
Fig. 4, the MRI converts to Alfven waves at k CT it [cf. Fig. 2(b)], and indeed the unstable 
Alfven waves show up for k > k CTit . The MRI itself appears to be largely unaffected by 
stratification, except that the maximum growth rate increases with increasing F z . It is 
unclear what would happen at higher fluxes because we are already close to the WKB limit 
([17]) at these high fluxes. Given that the vertical (&r = 0) MRI modes are unaffected by 
flux, they are not always the fastest growing modes in the presence of nonzero flux. Modes 
with a slight nonzero k^ in fact grow faster for the equilibrium parameters shown in Fig. 
4(c). The fact that F z can affect the MRI growth rates for kn ^ 0, albeit modestly, suggests 
that geometrically thick accretion flows with significant radial flux in the equilibrium could 
have a more marked effect on the MRI than the geometrically thin configurations considered 
in this paper. 

The photon bubble mode is always the most rapidly growing of the three MHD wave 
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instabilities that exist in stratified media. It appears to survive for k < k ciit , albeit at a 
smaller growth rate than expected from extrapolation of the high k behavior. Rotation 
therefore acts to mildly suppress, and does not kill, this instability when = 0. 

4.3. Nonzero Azimuthal Field 

We finally turn to the full dispersion relation with nonzero azimuthal field, but 
continuing to neglect gas pressure and vertical entropy gradients. Figure 5 presents results 
for the low wavenumber instability growth rates corresponding the the B^ = case shown 
in Fig. 4(b). As we stated above, the MRI is a robust instability which even in the presence 
of a nonzero azimuthal field and radiation, still exists for k < fc cr it. However, as expected 
from our results in section 3, the nonzero azimuthal field produces a marked reduction in 
the overall growth rate. The maximum MRI growth rates in Fig. 5 are comparable to the 
kji = growth rates shown in Fig. 1 at the chosen diffusion wavenumber k^\sv^ z /Vt = 0.1. 
The vertical radiative flux makes little difference to the MRI compared to the nonstratified 
case, just as we found in the previous subsection. 

The azimuthal field component has a noticeable effect on the other instabilities. First, 
the fast and Alfven waves are quickly stabilized. This is consistent with our findings in 
section 4.1, where we noted that azimuthal fields would stabilize the fast waves. For the 
parameters chosen in Fig. 5, equation fl55|) predicts fast wave stability for B^ > 0AB z . 

The photon bubble growth rate also decreases with increasing azimuthal field, in 
quantitative agreement with equation (Pf). At high values of B^/ B z , this mode appears to 
split into a high wavenumber part (the true photon bubble), and a low wavenumber part [cf. 
Fig. 5(b)]. The latter dominates the MRI growth rate except at the lowest wavenumbers, 
at least for the parameters chosen in Fig. 5. At high wavenumbers, this mode evolves into 
the Alfven mode. 



5. Astrophysical Relevance 

It turns out that diffusion is actually very important in standard a-viscosity models of 
accretion disks around black holes and neutron stars (Shakura & Sunyaev 1973). Neglecting 
relativistic correction factors, one can show that in the radiation pressure dominated inner 
zone, 

kdisH ~ 2a" 1 / 2 radians. (61) 
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Surprisingly, this relation is independent of accretion rate, radius, and central object mass. 
Combining this with the WKB condition for the MRI, k ~ Vt/vp^z ^ H" 1 , we find 



Depending on the assumed strength of the initial magnetic field, realistic accretion disk 
models (a ~ 0.01 — 0.1) will often be in the regime where photon diffusion can significantly 
slow the growth of the MRI in the presence of azimuthal field components (cf. Fig. 1). Our 
analysis shows that in this case the growth rate would be reduced to ~ Q(c g /vA<f>)- In the 
radiation pressure dominated inner zone, one can conceive of ^ being somewhat less than 
c r . Neglecting relativistic correction factors once again, 



where M is the mass of the central object. The MRI growth rate can therefore become 
quite sluggish if the azimuthal field energy density exceeds the thermal energy density in 
the gas by an amount comparable to the radiation/gas thermal energy density ratio. We 
stress however that even in this regime, the MRI still exists as a dynamical instability, 
and will still presumably grow eventually to nonlinear amplitude. It might be interesting 
to see if the resulting turbulence depends on how the initial magnetic field energy density 
compares to the gas and radiation thermal energy densities, but it is difficult to predict any 
results from linear theory alone. 

Note that all these conclusions refer to a fictitious, unstable initial state which cannot 
exist in nature. The equations above, in particular that for the diffusion wavenumber fl62|), 
depend on a, which is presumably a consequence of the MRI. It is therefore difficult to 
determine how the MRI will be affected when it is responsible for the ambient conditions in 
the disk in the first place. 

This fact is even more important when we consider the stratified calculations that we 
presented in section 4. Multiple instabilities are simultaneously present in this case, and 
the modes with the fastest growth rates presumably dominate the nonlinear evolution. One 
would at first conclude that the photon bubble/slow mode at high wavenumbers therefore 
controls the dynamics. From equation (0), the growth rate is at most 




(63) 




(64) 



substantially faster than the orbital frequency £1. 
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However, it is important to recognize the different sources of free energy which 
are driving these instabilities. The MRI feeds off the differential rotation in the flow, 
transporting angular momentum outward which ultimately allows material to accrete 
and release gravitational potential energy. We therefore view this instability as being 
fundamental to the dynamics of the accretion disk. Dissipation of the gravitational energy 
released and radiative transport processes set up vertical gradients which drive the photon 
bubble and other MHD wave instabilities. Convection will also presumably occur if the 
resulting entropy gradient has the right sign (Pietrini & Krolik 2000). It is the vertical 
gradients which provide the source of free energy for these extra instabilities, but these 
gradients are fundamentally determined by the action of the MRI. 

We speculate that on the largest vertical scales, where diffusion is negligible on an 
orbital time, the MRI will always operate as usual, launching perhaps a turbulent cascade 
to smaller scales. How this energy is ultimately dissipated is still unclear. If there is 
sufficient vertical flux, it may not be as simple as radiative damping of compressive MHD 
waves (Agol & Krolik 1998), as those waves are themselves unstable for some propagation 
directions. Given the high growth rates at short wavelengths from equation above, it 



seems likely to us that photon bubbles will play an important role in thermal transport on 
small scales. 

Concrete answers to the question of how radiative diffusion and radiative flux really 
affect the turbulent angular momentum and heat transport in the disk will require 
investigation by numerical simulations. One must be careful to choose the fictitious unstable 
initial state with care when pursuing such simulations, keeping in mind that accretion 
power, and therefore the MRI, is the driver for all this complex dynamics. 



6. Conclusions 

We may summarize our conclusions as follows. In the absence of an azimuthal field, 
the vertical (ka = 0) MRI modes have the standard, incompressible behavior elucidated 
by Balbus & Hawley (1991). If the azimuthal field component dominates the vertical field 
component and the Alfven speed is larger than the effective sound speed [cf. eq. (@)], then 
the growth rate of the MRI is reduced to ~ Q(c g /vA^>), well below the usual rapid growth 
rate of order the orbital frequency Q. However, the instability still persists. Large vertical 
photon fluxes in the equilibrium do not significantly affect the most rapidly growing MRI 
modes. 

At short wavelengths where rotation is unimportant, both the fast and slow MHD 
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waves are destabilized by sufficiently high vertical photon flux in the equilibrium. These 
instabilities require that the W£LV6S prop£lgclt6 clt some nonzero angle to the vertical, and 
affects upward propagating fast waves and downward propagating slow waves. The slow 
waves appear to always have the larger growth rates. At smaller wavenumbers these are 
well- described by Gammie's (1998) photon bubble dispersion relation fl49|) , but at larger 
wavenumbers they asymptote to a maximum value given by equation (|54j) . At small 
wavenumbers where rotation is important (k < k CTit ), the photon bubble growth rate is 
somewhat suppressed by rotation, at least at modest fluxes. Nonzero azimuthal fields 
stabilize the fast wave instability and reduce the photon bubble growth rate. In addition to 
the fast and slow wave instabilities, the vertical radiative flux also drives unstable modes at 
intermediate wavenumber which are related to Alfven waves at high wavenumber. 

The ambient conditions in the radiation pressure dominated inner region of standard 
alpha disk models are generally in a regime where radiative diffusion can slow the growth 
of the MRI. In addition, photon bubble modes are expected at short wavelengths with 
growth rates that can significantly exceed the orbital frequency. It is not clear how these 
instabilities will evolve into the nonlinear, turbulent regime, but we speculate that the MRI 
will dominate the largest scales and that photon bubbles (and perhaps unstable fast and 
Alfven waves) may affect both the dissipation and thermal transport on smaller scales. A 
complete understanding of the turbulent angular momentum and heat transport will require 
numerical simulations. 

Finally, we stress that we have concentrated here on geometrically thin accretion 
disk configurations where the photon flux is primarily in the vertical direction and the 
rotation curve is near-Keplerian. Geometrically thick or quasi-spherical, radiation pressure 
supported flows will have radiation pressure gradient components in the radial direction 
and may therefore exhibit more interesting dynamical coupling between the radiative flux 
and the MRI. 

We acknowledge useful conversations with Eric Agol, Yuri Levin, and Jim Stone. This 
work was supported by NSF grants AST-9970827 and PHY-9907949, and NASA grant 
NAG5-7075. 



A. Appendix A: The Zero Frequency Mode of the General Dispersion 

Relation 

Despite the fact that the perturbation equations (124|)-(|34|) and (|36| ) have eight separate 
tu-dependent terms, we only end up with a dispersion relation that is seventh order. The 
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reason is that there is a zero frequency mode, which we discuss in some detail in this 
appendix. 



Setting 10 = in the perturbation equations, the radial induction equation ( p3|) 
immediately implies that 5v R must vanish. The continuity equation (|24) and azimuthal 
momentum equation (|26|) then imply that 5t> 2 and must also vanish. The gas energy 
equation (|28|) is then trivially satisfied, which is the mathematical reason why a zero 
frequency mode is a possible solution of the perturbation equations. 

Excluding the trivial azimuthal flux equation ([IT]), the remaining seven equations may 
be written 

- 2np6v4 + ik R Up + ^6E^j + i^B z 5B z - ^BJB R = 0, (Al) 

ik z (6p + -6e) + g5 P = 0, (A2) 
k • 5F = 0, (A3) 
5F R + k R 5E = 0, (A4) 



8F z + F z ^ + -^k z 5E = 1 (A5) 
- R^z5B R - ih z B z bvs = 0, (A6) 



and 

K5B R + k z 5B z = 0. (A7) 

These seven equations have eight unknowns, so that any one of the nonzero perturbation 
amplitudes determine the other seven for this zero frequency mode. 

It is perhaps surprising that this mode has such a complicated eigenf unction, and it 
is therefore worthwhile understanding its physical meaning. Equation ( |A~6|) is particularly 
interesting. This arose from the azimuthal component of the induction equation which 
describes, for example, how a radial magnetic field component is sheared out by the 
differential rotation to produce a time-dependent azimuthal field component. In this case, 
however, this is exactly compensated by an azimuthal velocity perturbation which creates 
azimuthal field in the opposite direction out of vertical field. 

Further insight into the physical nature of this mode may be obtained by taking various 
limits. In the limit k z — > 0, we have 5B R = 5p = 5E = 8F = 0, and the only nontrivial 
condition comes from the radial momentum equation: 

- 2QpSvA + ik R 5p + i^-BJB z = 0. (A8) 

4:71 
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In other words, any random field of perturbations of Sv^ can be balanced by radial gradients 
of gas and/or magnetic pressure to produce a slightly different equilibrium. 

In the opposite limit k,R — > 0, we have 5v<f, = 5B Z = 5Br = 5F Z = 0, and 

ik g (Sp+^SE^+gSp = 0, (A9) 

5F R + - k R 5E = 0, (A10) 

F,^ + ^-W£ = 0. (All) 

In this limit, a random field of density perturbations is balanced by vertical gradients in gas 
and radiation pressure, again producing a slightly different equilibrium. 

Finally, we note that if we neglect the vertical stratification of the equilibrium, then 
this zero mode has all perturbations vanishing except for the density. In this case a random 
field of density perturbations will merely be advected along with the equilibrium flow 
without producing additional net accelerations. 
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Fig. 1. — Maximum growth rate of the MRI in a medium with no stratification, as a function 
of azimuthal field strength and diffusion wavenumber. We have adopted a Keplerian rotation 
curve (k = Q) for these calculations, and have set c g = 0, c t /va z — 10, and = 0. Each 
curve corresponds to different values of v^/v^z, which are labelled. Note the substantial 
reduction in growth rate at modest azimuthal field strengths and low diffusion wavenumber. 
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Fig. 2. — Real part of the mode frequencies in a medium with no stratification and a vertical 
magnetic field (B^ = 0), as a function of total wavenumber. We have assumed a Keplerian 
rotation curve and k z /k = 0.8, c t /vaz = 10, and c g /vAz = 0.2. Each figure corresponds 
to a different value of kdisVAz/^'- (a) 10, (b) 0.1, and (c) 0.01. The critical wavenumber 
fc cr it = [— R{d£t 2 / dR)] 1 / 2 / vaz above which the MRI is stabilized is indicated in each figure. 
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Fig. 3. — Imaginary part of the mode frequencies in a stratified medium with vertical 
magnetic field (B^ = 0) as a function of total wavenumber. The rotation curve is Keplerian 
and k z /k = 0.8, c r /v\ z = 10, c g /vj^ z = 0.4, k A ^v^ z /K = 0.1, and F z /(3pv\ z ) = 1000. The 
dashed curves show the analytic approximations (|49|), (|53|), and (0). [We have neglected 
the Brunt- Vaisala terms in the dispersion relation (|4"6D and also set g/(nvAz) = 10 and 
(g/ K 2 )dln p/dz = —1. However, these last two parameters do not significantly affect the 
results and could just as well have been set to zero. This is consistent with the WKB 
approximation.] 
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Fig. 4. — Imaginary part of the mode frequencies in a stratified medium with vertical 
magnetic field (B^ = 0) as a function of total wavenumber. The rotation curve is Keplerian 
and k z /k = 0.8, c r /v Az = 10, k diS v Az /K = 0.1, and F z /{3pv% z ) = (a) 30, (b) 200, and (c) 
1000. Gas pressure has been neglected and the vertical entropy gradient has been set to zero 
(i.e. N T = 0). 
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